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A COMPARISON OF TECHNIQUES 
FOR INVERSION OF RADIO-RAY PHASE DATA 
IN PRESENCE OF RAY BENDING 

By H. Andrew Wallio and Mario D. Grossi* 

Langley Research Center 

SUMMARY 

Derivations are presented of the straight-line Abel transform and the seismological 
Herglotz-Wiechert transform (which takes into account ray bending) that are used in the 
reconstruction of refractivity profiles from radio-wave phase data. Profile inversion 
utilizing these approaches, performed in computer -simulated experiments, are compared 
for cases of positive, zero, and negative ray bending. For thin atmospheres and iono- 
spheres, such as the Martian atmosphere and ionosphere, radio-wave signals are shown 
to be inverted accurately with both methods, For dense media, such as the solar corona 
or the lower Venus atmosphere, the refractivities recovered by the seismological 
Herglotz-Wiechert transform provide a significant improvement compared with the 
straight-line Abel transform. 


INTRODUCTION 

Radio occultation measurements of planetary atmospheres and ionospheres have 
been an integral part of the standard scientific investigations performed with planetary 
space probes since the Mariner series of the mid-1960's. The Mariner IV spacecraft 
provided the first reliable measurements of the atmosphere and ionosphere of Mars by 
this technique in 1965. (See refs. 1 to 3.) Similarly, the Mariner V spacecraft provided 
refractivity profiles for Venus. (See refs. 4 to 6.) 

One of the basic observables used in these experiments is the observed frequency 
of the radio signal received from the spacecraft. This frequency is compared with a 
fixed stable oscillator to derive an observed Doppler signal. By subtracting the amount 
of shift expected from the spacecraft-Earth link geometric changes from the observed 
Doppler signal, a "Doppler shift residual" is determined. For situations in which both 
link terminals are outside the medium being probed, the presence of a nonzero residual 
indicates a change of the emitted and/or received angle of the ray that propagates between 
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the links and through the atmosphere of the planet being probed. The inversion problem 
is to obtain the atmospheric refractivity profile as a function of radial distance from 
the planetary center of mass from the observed Doppler shift and/or Doppler shift 
residuals. From the refractivity profile and certain assumptions about the constituents 
of the atmosphere and the planetary magnetic field, such atmospheric properties as the 

electron number density and its distributions, pressure, temperature, density, neutral 

/ 

number density, and their radial dependence can be deduced. 

In the past, the inversion of the data has been performed with a variety of geometric 
optics techniques ranging from the closed-form Abel transform (based on the assumptions 
that the medium possesses spherical symmetry and based on the additional hypothesis-, 
generally accepted by occultation experimenters, that there is straight-line ray travel 
(ref. 7)) to ray trace model-fitting approaches (based on iterative model update schemes 
(ref. 8)). For a planet with a thin atmosphere, such as Mars, all these methods are, in 
principle, adequate approximations to the ray-propagation theory. For dense media, like 
the atmospheres of Jupiter and Venus or the Sun corona where ray bending is large, 
closed-form inversion algorithms which utilize a straight-line ray approximation are 
inadequate. A suitable method to solve the problem in these cases is available from 
seismology and is based on the Herglotz-Wiechert approach for interpreting seismic 
data. (See refs. 9 and 10.) 

In seismology, this bent-ray-path technique leads to the derivation of velocity- 
depth profiles and to the calculation of the position and the identification of geological 
features in the Earth's mantle. The observables are the travel times of seismic waves 
between stations located at known positions on the Earth's surface. Phinney and 
Anderson (ref. 11) have shown that this method is applicable to situations where the 
observables are the Doppler residuals of the radio occultation measurements. This 
method has been independently derived, programed, and applied to the Martian atmo- 
sphere (ref. 12) and this report is an extension of this earlier work to include denser 
atmospheres. 

The purposes of this paper are to present the equations of the closed-form straight- 
line Abel transform and the Herglotz-Wiechert transform, and to compare the refrac- 
tivities and minimum probing radius recovered by the two transforms for various cases 
of ray bending. 


SYMBOLS 

A amplitude factor of electric field wave, V-m"* 

a,x,Z limits of integration 
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c speed of light, km-sec"-*- 

D Doppler shift, sec - * 

*\tm Doppler shift due to atmosphere, sec - * 

D g Doppler shift of unperturbed radio ray, sec - * 

dl differential actual radio ray-path length, km 

dl' differential straight-line radio ray-path length, km 

dr differential of r, km 

dS differential of S, km 

d0 differential of 0, deg 

E electric field vector, V-m'* 

f frequency of radio wave, sec'l 

G phase factor of electric field wave 

imaginary part of complex G 
G r real part of complex G 

g(x), g’(x), £, ii dummy functions used to illustrate Abel's integral equation 
i incidence angle, deg 

K propagation constant, km - * 

K 0 propagation constant in free space, km - * 

L total path length, km 


l 


actual radio ray-path length, km 



straight-line radio ray-path length, km 
refractivity, 10®(n - 1) 


refractivity at minimum radial distance of radio ray 


index of refraction 


index of refraction at minimum radial distance of radio ray 


impact parameter, km 


radial distance from planetary center of mass, km 


starting position of radio ray, km 


radial distance of minimum point of radio ray, km 


position of spacecraft on trajectory, km 


spacecraft velocity along trajectory projection, km -sec' 


T,t,x defined dummy variables 


„ _ dT(p) 


geometric angle (see fig. 3), deg 


variable between 0 and 1 used in Abel’s integral equation 


deviation angle, deg 


inductive capacity, N"*(C-m"l)^ 


defined variable, km 


starting value of 77 for a given ray, km 


minimum value of 77 for a given ray, km 
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6 unit of angular measure (see fig. 2), deg 

X wavelength, km ■* 

p perpendicular miss distance, km 

<p angular measure of phase cycles 

<£ a defined ray residual, cycles 

0 U geometric phase term, cycles 

\f/ radio-wave emission angle (see fig. 3), deg 

i// unperturbed radio-ray emission angle (see fig. 3), deg 

Subscripts: 

0 starting position 

L last position of integration 


ANALYSIS 

Straight -Ray Abel Inversion Transform 

The differential phase -path length is the difference between the straight-line geo- 
metrical distance between transmitter and receiver and the phase-path length for radio 
waves and is obtained from an integration in time of the Doppler residuals. The integra- 
tion constant is zero when provisions are made for starting the integration from a posi- 
tion time where the probing link is completely external to the medium under investigation. 

Figure 1 depicts the straight-ray Abel transform geometry. By assuming that the 
index of refraction is a spherically symmetric function, 

pOC 

Differential phase-path length = \ dl’ - \ n(r) dl (1) 

- 00 J - 00 

By assuming that the radio-wave ray path is a straight line and sufficiently close to 
the geometric straight line that ctf ~ dZ', 
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c 00 

Differential phase-path length = \ fl - n(r)l 

- 00 


dZ' 


( 2 ) 


The integral is a function of the distance of the ray from the center of the planet at 
the point of its closest approach. 

Define 


T (p)= f [l - n(rj] dZ’ 

- OO 

Then, by changing the variable dZ' (see fig. 1), 

_ r”i»Wr_dr 

J r 2 - p2 .La I? 


(3) 


T(p) 


\jr 2 - p2 


(4) 


In this equation the radicand in the denominator is nonnegative, because for every value 
of the miss-distance range p of the ray from the center of the planet, r > p. 

Let t s and 2r dr = -4r dt. Equation (4) becomes 


o 10-6 V dt 

T(p) = j XVt ' 


1/ P 2 


t2 /r 




, 9 10” 6 N^— \ dt 

i/p 2 [ft] 




(5) 


Let x = and p = w— ; then equation (5) becomes 

O" i x 


10 _6 Nf-M d t 

■x l/t. 


JL T /JJ\ - f - W 
>/5 w Jo t Wf~r t 


(6) 


Applying Abel's integral equation in the form given in reference 13 as if 


g(x) = J 


,x m(I) dg 
a (x - %)P 
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then 


_ sin f3n C Z g'(x) dx 

77 4 (Z - x) 1 '^ 


(under the appropriate conditions as given in ref. 13) to equation (6) results in 


\]ftj _ _ 106 r»t dx jj/x 

71 *^0 \jt - x 


Let t s — and then equation (7) becomes 

r*" 


t/-M 


- r 3 N(r) = io! r 1/f2 dx lg 
77 4 /I 


aq_ T mi= -i 

d* _/x \\/x/_ 2 


Tf-1) T’(4 

\|fx/ 1 \/x 

x 3/2 2 2 


so that equation (8) becomes 


W w 

.3,^. ID 6 

2^ Jo n 


where T’ = Let x = ~ and 2p dp = dx; then equation (10) becomes 


*■» ■ -£ II 


+ T’(p) 2p dp 


I r 2 p2 
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or rewritten 


N(r) = 10« r PT(P) + P 2 T'(p) ^ 

70*2 \jp2 _ r 2 


(12) 


and 


, 1 f°° pT(p) + P 2 T’(p) 

n(r) = 1 + — — === 

nr* J r \jp^ . 


dp 


(13) 


Because of the change in the limits of integration, the radicand in the denominator is 
nonnegative. In fact, for every radial height r, N(r) and n(r) are obtained from 
columnar measurements made at miss distances p always larger than r. 

Thus, under certain restrictive assumptions, the index of refraction of the medium 
being probed is shown to be a function of the differential phase-path length and its deriva- 
tive which are known observables. 


Herglotz-Wiechert or Seismic Inversion Transform 

In seismology, the observable is the central angle 6 subtended by seismic wave 
rays in the Earth. (See fig. 2.) The analytical steps are arranged in such a way that the 
velocity-depth profiles are obtained by operations on 6. 

In radio occultation measurements, the observables are the Doppler shift residuals 
so that the theory used in seismology must be modified so that the refractivity or index - 
of-refraction profiles are obtained by operations on the Doppler shift residuals. The 
important aspect of the Herglotz-Wiechert approach is that the straight-line ray approxi- 
mation is not required. 

The geometry of the occultation is shown in figure 2. The index of refraction n is 
assumed to be a radially dependent function, and from Fermat's principle, the first vari- 
ation of the phase-path length should be zero, 


6 


I 


n dl = 0 


(14) 


where d = dr^ + dd The Euler-Lagrange equation which yields a minimum for 
equation (14) is 


nr 


2 M 
dr 


1 + r 2 y 
L Vdr/ 


U/2 


= Constant 


(15) 
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and from the boundary conditions at the minimum radius, Snell's law in spherical geom- 
etry is determined, that is, 


p = nr sin i = nprp h Impact parameter 

The impact parameter is a constant for a given ray. 
Define the variable t) as 


r\ = nr = — — r 
1 sin 1 


From the definition of the path length 


r di = r 

[i + r2(^) 2 ] 

j j 

[_ \dr/ J 


1/2 


dr 


(16) 


(17) 


(18) 


and the application of equations (15) to (17), the integral function for the path length is 
found to be 


L = 2 



(19) 


The phase (angular measure) along the ray path is 


= | I n dZ 

and by using equations (17) and (19), 

<Mp) = 2 | J ° -f! 7 ? 2 - p 2 ) dr 

r P 

By applying the operator L to equation (21) where 



(20) 


( 21 ) 


( 22 ) 


and where t?q = = rjp by integrating over p, by interchanging the order of integration, 

by integrating by parts repeatedly, and by performing a final integration over the appro- 
priate regions, the phase function (eq. (21)) can be shown to be inverted (see appendix) to 
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(23) 



Thus, the refractivity is found as 



where <p u is a geometric term, then 



where x is a dummy variable such as time to occultation or satellite position. 

Thus, if the impact parameter p (eq. (16)) can be determined, the minimum 
probing radius and the refractivity at that radius are determined for each radio-wave ray 
as a function of the known differential phase path by equations (23) and (26), respectively. 


Determination of Impact Parameter 

In order to evaluate the impact parameter, a ray-optical treatment is considered 
from an Earth-centered coordinate system. The spacecraft is thought of as sampling 
phase along an arc through a family of constant phase surfaces and the radio waves are 
emitted orthogonally to these surfaces. 

Then the equation for the directional derivative (ref. 14) is 

Sir . vd> = ^ (27) 

dS r dS 

where dr is the tangent line to the projection of the spacecraft trajectory in the plane 
of the spacecraft, planet, and Earth at any point S along the trajectory, dr/dS is a 
unit vector, and V<£ is the direction of propagation. 
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Therefore, the angle between the radio ray and trajectory becomes (see fig. 3) 

cos^ = ^|v<^j _1 (28) 

In order to evaluate the radio-ray emission angle from the directional deriva- 
tive, the gradient of the phase V0 or the eiconal equation must be determined. To 
determine this equation, assume a nonconducting isotropic medium where the inductive 
capacity e is a function of position. Then the wave equation for the electric field vector 
becomes 

V 2 E + K 2 E = -V^E • = V 2 E + Ko 2 n 2 E (29) 

where Kq = f/c is the propagation constant measured in cycles. If it is assumed that 
the spatial change of e is small compared with a wavelength 


Ve 

e 



(30) 


then the wave equation reduces to its homogeneous form 
V 2 E + K 0 2 n 2 E = 0 


Assume that the form of the solution is 


E 


= Ae- iK ° G 


(31) 


(32) 


where G = G r + iGj so that surfaces of G r = Constant are surfaces of constant phase 
and surfaces of Gj = Constant, are surfaces of constant amplitude. Substitution of equa- 
tion (32) into equation (31) yields 


V 2 A + K^aJti 2 - (VG) 2 J - iK Q [AV 2 G + 2 (va) • (vg) 
If Kq is large in the sense that 


= 0 


(33) 


V 2 A « Kq 2 


(34) 


and 


AV 2 G + 2(va) * (vg) « K 0 


(35) 
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(therefore regions of diffraction, focal points, caustics, and sources are excluded), then 


(VG) 2 = n 2 


(36) 


Now K 0 G = (p + ft for a real r? so that 

(V0) 2 = n) 2 (37) 

which is the eiconal equation. 

Since n is a function of position if the spacecraft shown in figure 3 is inside the 
planetary atmosphere, the profile of the atmosphere must be known before equation (37) 
can be solved for | V0| . 

If, however, the spacecraft is above the atmosphere, then n becomes unity and 
equation (37) becomes 

(V0) 2 = (|) 2 (38) 

Substituting equation (38) into equation (37) yields 

cos ^ = f (39) 


By letting time be the independent variable and D signify an instantaneous 
Doppler shift, 


d(ft _ i dcp _ p 

dS s * S 

Then D can be thought of as being composed of two parts 


(40) 


D - D e + D a j. m (41) 

where D e is the expected geometric Doppler shift from the predicted trajectory and 
D a tm is the residual normally attributed to the atmosphere. Therefore, 


cos 'P = ^ ( D e + D atm) 


(42) 


If there is no atmosphere 


cos <// p = D 
e fS 1 


(43) 
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where \p e is the angle between the trajectory tangent and the unperturbed ray or 
straight-line distance between the spacecraft and Earth. 

Algebraic manipulation of equations (42) and (43) yields the deviation angle 6 
(fig- 3) 


or 


5 = p - \p e 



°atm) 2 ] 1/2 - (D s 2 - D e 2 ) 


1 / 2 ' 


D atm + 2D e 


6 = 2 tan - * 




D + D e 



( 44 ) 


(45) 


where 

D observed Doppler shift 

Datm Doppler residual due to atmosphere 

D e Doppler shift of unperturbed ray 



Thus, the impact parameter for a ray is 

p = r s sin (a e - 6) (46) 

where 

r s distance of spacecraft from planet center of mass 

of e a known geometric angle (fig. 3) 

Implied in equation (46) is the dependency of each ray impact parameter on either of the 
related variables of spacecraft position or time. 
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RESULTS AND DISCUSSION 


A comparison is made of the recovered refractivity profiles obtained by using the 
straight-line Abel and Herglotz-Wiechert transforms. Three types of atmospheres have 
been used to compare the two forms of transforms: Venus (where the radio rays are bent 
positively inward by the neutral atmosphere); Mars (where ray bending is negligible); and 
the solar corona (where the rays are bent negatively outward by the electron plasma). 

For this study it was assumed that the Earth-planet distance was sufficiently large 
to insure that the radio rays emerging from the atmosphere can be considered parallel. 
Therefore, all the observable atmospheric frequency shifts can be considered to arise 
from effects that have occurred between the spacecraft and an imaginary reference plane 
located between the Earth and the planet and beyond the sensible atmosphere. The phase 
data for the three models were generated, for a frequency of 2000 MHz, by a ray-tracing 
program; this program was a modified version of a program reported in reference 15. 

The assumed Venus atmospheric model has a relatively small ionosphere component 
compared with the neutral layer. For this model, ray bending in the lower atmosphere is 
significant and the differences in the recovered parameters between the straight-line Abel 
and the seismic Herglotz-Wiechert approach are clearly noticeable. The recovered 
refractivity from the straight-line Abel (Abel), the seismic Herglotz-Wiechert (Seismic), 
and the model refractivity (Model) against the appropriate planetocentric radial position 
from Venus are plotted in figure 4 for the measured phase function. 

In table I is a listing of the differences in the recovered refractivity and minimum 
probing radii of the two transform methods from the respective actual refractivity and 
minimum probing radius. From table I and figure 4, it can be seen that as the radio 
rays probe lower into the atmosphere and ray bending increases, the divergence of the 
recovered values increases, the seismic Herglotz-Wiechert method giving consistently 
better recovered values. 

The Martian atmosphere model is similar to the one developed from the Mariner IV 
results. For this example, ray bending is small and the straight-line Abel transform has 
been shown in the past to be a fully adequate approximation. The recovered refractivity 
from the straight-line Abel transform (Abel), the seismic Herglotz-Wiechert transform 
(Seismic), and the model refractivity (Model) against geocentric radius of Mars are 
plotted in figure 5 for the measured phase function. 

Table II is a listing of the differences in recovered refractivity and minimum 
probing radius of the two transform methods from the respective actual refractivity and 
minimum probing radius for Mars. From figure 5 or table H, it can be seen that the 
errors in the recovered refractivity of the two transforms are about the same. Also, the 
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minimum probing radius found by each transform is very close to the actual minimum 
probing radius because of the small amount of ray bending. 

A Baumbach electron density model was used for the solar corona. For this case 
ray bending is again very large and the deviations between the classical straight-line 
Abel transform and the seismic Herglotz-Wiechert transform are clearly noticeable. 

The recovered refractivity from the straight-line Abel (Abel), the seismic Herglotz- 
Wiechert (Seismic), and the solar model refractivity (Model) against their respective 
heliocentric radial position are plotted in figure 6. 

Table HI is a listing of the differences in recovered refractivity and minimum 
probing radius of the two transform methods from the respective actual refractivity and 
minimum probing radius for the sun. From figure 6 or table HI, it can be seen that the 
divergence of the recovered refractivities and minimum probing radii compared with the 
actual model values increases as the ray bending increases, the seismic Herglotz- 
Wiechert transform yielding the parameter values closest to the model values. 

CONCLUDING REMARKS 

The work reported herein has shown that profile inversion of phase data collected 
from strongly bent radio rays requires the use of techniques like the seismological 
Herglotz-Wiechert transform which do not rely on the validity of the straight-line 
approximation. 

For highly refractive media like the solar corona or the lower Venus atmosphere, 
the recovered refractivities by the Herglotz-Wiechert transform provide significant 
improvement compared with the straight-line Abel transform (a reduction of refractivity 
error at the lowest probed altitude from over 60 percent to less than 20 percent for the 
solar corona and a reduction of error from 150 percent to 3 percent for the Venus lower 
atmosphere). 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., October 20, 1972. 
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APPENDIX 


TRANSFORMATION OF PHASE EQUATION TO EQUATION 
FOR MINIMUM PROBING RADIUS 


This appendix gives an outline of the mathematical steps used to transform the 
phase function (eq. (21)) to the equation for the value of the refractivity at the minimum 
point of the probing ray. 

If the phase function (from eq; (21)) 

«p)=2|^°3!(, 2 -p2)- 1/2 | d , (Al) 

is given, let rj be such that Tq = r^ = R, then 7 ?q = ^ Apply the kernel 

\ _ 1 /2 

p(p^ " Vi) to equation (Al) and integrate over p from rjj to ?]q to yield 



P0 




dp = 2 


f 

c 




(A2) 


Interchanging the order of integration results in 


r% 

P (p dp f 

p's!**, 

r 7 ! 0 

2p dp 

^1 

lt 2 -V c 

\ r *> v 


/(p 2 - vi 2 ){v 2 ~ P 2 )_ 


f C 71 0 o dlog e r , 

-o) ni —&r~ " 


which reduces to 

lQ g e r _ c d /r 770 p<a ^ \ (A 

dT] rj^irf dri {j^ ^2,^2 J 

Let ri — and integrate both sides of equation (A4) over p. from t?q to 77 to yield 
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APPENDIX - Continued 



Let ri^ = ji and then integrate equation (A5) by parts to yield 



Interchanging the order of integration and integrating with respect to ju result in 



which reduces to 



0^p2 - rft $ 
V 2 P P 2 



(A8) 


The second term on the right-hand side of equation (A8) becomes, after integrating 
by parts, 



The third term on the right-hand side of equation (A8) becomes, after integrating by 
parts, 



Combining equations (A10) and (A9) in equation (A8) yields 
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APPENDIX — Continued 


log f 


~Hv) 

j 

C Vo f 0P 4>P 

1 

Pnc — ■ 

l (n) 

L r o. 

rf L 

V \p 2 \] p 2 - rft 7|2^ p 2 _ ^2 

‘ V 


W 


1# 

P 


U cosh " 1 (f) ■ I C 0 S 'Kp)w 



Equation (All) reduces to 


r(v) = r 0 exp 


3 r c ° sh " i ®5 * dp 


^0 


and since r(rj) = rg, the negative value of double- valued cosh'^^j must be taken. 
At the "turning point" in the atmosphere (where sin i = 1) 


rjp = p = r 0 sin i Q = n p r p 

Therefore, transforming equation (A12) to a dummy variable dependence and using 
equation (A13) 



cosh 


-1 


P(x) 


P( x l) 


1 ^ dx' 

p[xj dx 


Using equation (A13) again yields 


n p = sin i 0 




cosh 


-1 


r p(x)i 

1 d (p 

_P( x l)_ 

p(x) dx 



Now <p = <p^+ <fi^ so that if there is no atmosphere, <p = cf> and equation (A15) 
becomes 


1 



cosh 


-1 


p(x) 

_p( x l) 


1 d< ^M 
p(x) dx 



and therefore 


(All) 


(A12) 


(A13) 


(A14) 


(A15) 


(A16) 
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APPENDIX - Concluded 


7 


N p = 10' 


exp 


V 



cosh" 


i r p(x) 

1 d ^a 

_P(x L )_ 

p(x) dx 



(A17) 


which is equation (26). 
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Figure 2.- Seismic transform geometry. 
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Figure 3.- Impact parameter geometry. 
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Figure 4.- Model of Venus refractivity profile and its reconstruction by Abel transform and 
Herglotz-Wiechert transform. (Rectangles traced around the pairs of circles, squares, 
and triangles identify the original and the reconstructed point.) 




Figure 5.- Model of Mars refractivity profile and its reconstruction by Abel transform and 
Herglotz-Wiechert transform. (Rectangles traced around the pairs of circles, squares, 
and triangles identify the original and the reconstructed point.) 
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Figure 6.- Model of Sun corona refractivity profile and its reconstruction by Abel transform and 
Herglotz-Wiechert transform. (Rectangles traced around the pairs of circles, squares, and 
triangles identify the original and the reconstructed point.) 
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